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Recently a generalized master equation was derived that extends the Lindblad theory to highly 
non-Markovian quantum processes (H.-P. Breuer, Phys. Rev. A 75, 022103 (2007)). We perform a 
stochastic unravelling of this master equation by considering n random state vectors that satisfy the 
corresponding stochastic differential equation for a piecewise deterministic process. As an applica- 
tion we consider a two-state system randomly coupled to an environment consisting of two energy 
bands with finite number of levels. Our numerical results are compared to results obtained from the 
time-convolutionless (TCL) projection operator method using correlated projection superoperators 
and the exact solution of the Schrodinger equation for this system. 

PACS numbers: 03.65.Yz, 42.50.Lc 



I. INTRODUCTION 

The success of present and future quantum technolo- 
gies relics almost entirely on the quantum device's in- 
teraction with the environment it is in. Decoherence 
and dissipation phenomena dictate how much informa- 
tion can be transmitted from one quantum manipulation 
to the next. Decoherence, which is the loss of phase co- 
herence between superpositions of quantum states, and 
dissipation, which is the leakage of population from the 
system to the environment, are major hurdles to the re- 
alization of realistic quantum technologies. As a result, 
the investigation of the dynamics of open quantum sys- 
tems, is of utmost importance to our understanding of 
such undesirable phenomena [l[. 

Most approaches to the investigation of open quan- 
tum systems are based on Markovian assumptions, which 
makes use of the Born and Markov approximations that 
ultimately lead to the quantum Markov equation in Lind- 
blad form In most cases this Lindblad master equa- 
tion is stochastically unravelled enabling the efficient use 
of stochastic wave function methods to analyze the dy- 
namics of the open quantum system. These methods 
have prominence in applications to many quantum opti- 
cal systems [1, S 0, H 111 • 

In some instances however, open quantum systems as- 
sociated with more realistic quantum technological pro- 
cess are classified as non-Markovian. Some prime in- 
dicators of the presence of non-Markovian effects and 
the failure of Markovian approximations are when the 
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system-environment couplings are strong or when the ini- 
tial states are classically correlated or entangled. Some 
examples of non-Markovian systems include spin star sys- 
tems H, and circuit QED [ll|, 111. Various tech- 
niques have been developed to describe non-Markovian 
quantum process. Generalized or non-Markovian mas- 
ter equations have first been introduced in Refs. [3, 
The Nakajima-Zwanzig formalism p^ . [l6| and the time- 
convolutionless projection operator method [13, [H, fioj 
have proved to be useful in deriving approximations 
based on projection operator techniques. The latter 
method, employing correlated projection superoperators, 
was recently used to derive a non-Markovian generaliza- 
tion of the Lindblad equation . Stochastic wave func- 
tion methods have also been proposed and developed for 
non-Markovian quantum master eq uations [2ll . [2^ . [23| 
and more recently by Piilo et al [24| . 

In this paper we perform a stochastic unravelling of the 
generalized Lindblad master equation which allows for 
the use of traditional Markovian stochastic wave- function 
simulations. This approach is applicable to both time 
dependent and time independent rates. As an applica- 
tion we consider a two-level system coupled to an en- 
vironment consisting of two energy bands, each with a 
large number of energy levels. Due to its highly non- 
Markovian characteristic, this model has gain ed some 
interest over the past couple of years pol. l25l. [26l . 
In Ref. [2^, the time-convolutionless projector opera- 
tor technique and the Hilbert-space-average method was 
used to analyze this model; our Monte Carlo simulations 
are compared to the former technique. Similar models of 
this type have also been studied before. These include 
the model by Esposito and Gaspard [11] , and the models 
by Bixon and Jortner [2^ in the late sixties [13] ■ 

Huang et al [Slj have recently discussed an unravelling 
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for the generalized Lindblad equation as applied to the 
model being discussed in this paper for the case of con- 
stant rates. Here, we are interested in the case of time 
dependent rates involved in the strong coupling regime 
of this model. 

The paper is organized as follows. In Sec. |TT] we de- 
scribe the stochastic unravelling of the generalized Lind- 
blad equation that was derived in Ref. [201. Sec. lIIII we 
describe the model used and quote results obtained from 
the TCL expansion using correlated projection superop- 
erators as derived in Ref. [1^. In Sec. |TV] we perform 
the stochastic wave-function simulations for the model 
and consider both the weak coupling and strong coupling 
regimes. Results and conclusions follow respectively in 
the last two sections. 



II. THE GENERALIZED LINDBLAD 
EQUATION AND ITS STOCHASTIC 
UNRAVELLING 

The general form of the non-Markovian master equa- 
tion, obtained from the application of correlated projec- 
tion superopcrators. derived in Ref. [20| is given by 



dt' 

(1) 

where i,j = 1,2, ... ,n with _ff* being arbitrary Hermi- 
tian operators and arbitrary system operators (with 
?i = 1). This master equation preserves the normaliza- 
tion and positivity of the density matrix, /Oi(i). Following 
the procedures discussed in Ref. [l[, the stochastic un- 
ravelling of this equation is obtained by taking n random 
state vectors that satisfy the stochastic differential 

equations for a piecewise deterministic process in Hilbert 
space: 



dNiit). (2) 



The unnormalized density matrices pi are then deter- 
mined by the expectation values 



p,(i)=E(|^,;(t))(V',;(i)|). 



(3) 



The second term on the right hand side of Eq. ^ 
contains the Poisson increments dNl (t) which satisfy. 



and 



where 



dNldNl, ^5,,, 5 J J, dNl 
E{dNl) = Mldt 



(4) 
(5) 
(6) 



The first term on the right hand side of Eq. ([2|) de- 
scribes the deterministic drift of the process given by 

and with this, the deterministic pieces of the process are 
described by the differential equation 



The jumps are given by 




(8) 



(9) 



which occur at the rate M-j,. It should be noted that all 
state vectors jump simultaneously. 

Using the Ito calculus [l|, Is^l for piecewise determinis- 
tic processes, one can demonstrate that the expectation 
values given by Eq. ([3]) satisfies the generalized Lindblad 
equation ([T]). The stochastic unravelling nicely illustrates 
the fact that the master equation preserves the positiv- 
ity of the Pi since an expectation value of the form Q 
automatically represents a positive matrix. 

A further remarkable property of the piecewise deter- 
ministic process is that the total normalization is strictly 
preserved under the time evolution: 



(10) 



This implies that the trace of the reduced density matrix 

i 

is strictly conserved (not only on average): 

trpsit) = ^trp,(<) = ^(V'z(i)IV'^W> = 1- (12) 

i i 

Moreover, the quantities (V'ilV'i) can vary only between 
and 1 and the norm of all components is bounded. This 
means that there is no exponential growth of the norm 
of the state vectors as in other Monte Carlo approaches 
to non-Markovian quantum dynamics. 



III. THE MODEL AND RESULTS FROM THE 
TIME-CONVOLUTIONLESS METHOD 



We consider the two-state system coupled to an en- 
vironment consisting of two energy bands, each with a 
finite number of evenly spaced levels. The total Hamil- 
tonian in the Schrodinger picture is given by [25| . 

H = iASa, + ^-^ni|ni)(ni| 
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FIG. 1: A two-state system, with level distance AiJ, cou- 
pled to an environment consisting of two energy bands, each 
with a finite number of evenly spaced levels A'^i and A^2- 5e 
is the width of the bands and V is the system-environment 
interaction potential. 



where, the system-environment interaction potential has 
the form 

V^(ni,7i2) = A ^ c(ni, n2)cr+|ni)(n2| + /i.e.. (14) 



Here, riijni labels the levels of the lower(A^i 
levels) /upper(A'2 levels) energy band and A gives 
the overall strength of the interaction. 8e is the width of 
the upper and lower energy bands and AE' is the level 
distance of the two-state system. The coupling constants 
c(n\^ni) are complex Gaussian random variables with 
zero mean and unit variance. 

We consider the initial state where only the lower 
band is occupied. For the weak coupling case where 
8tt ^ 1, the second order of the TCL expansion using 
correlated projection supcropcrators, which we call new 
TCL2, gives the following equations of motion p5| : 

^Pl = 7iCr+p2CT- - Y{cr+cr-,Pl} (15) 
= 72Cr_pi(T+ - y{cr_CT+,p2}, (16) 

with the following solution for the population of the up- 
per level. 



P\\ = Pii(O) 
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71 



-(71+72)* 



71 +72 71 + 72 



(17) 



For the case where the times t do not satisfy the con- 
dition Set ^ 1 (strong coupling), the second order of the 
TCL expansion using correlated projection supcropcra- 
tors, which we call new TCL2(t), the equations of motion 
are: 



d 
di 



dt 



pi = / dtih{t - ti)[2jia+p2a^ 
Jo 

-72{<T+(T-,pi}], 



P2 



dtih{t - ti)[2720-_pi(T+ 
-7l{cr_cr+,P2}], 



(18) 



(19) 



where 72^1(^—^1) is the two-point environment correlation 
function such that 



h{t) 



Se sm^{5et/2) 
2^ (fct/2)2 ■ 



(20) 



The solution for the populations of the upper level in this 
case is given by. 



Pii = Pii(O) 



71 



71 +72 71 + 72 



(21) 



where 



r(t) = 2(71+72) / dh f dt2h{h-t2). (22) 



For both cases, the relaxation rates are given by, 

.2 



7l,2 



27rA2iVi.2 



Se 



(23) 



IV. 



STOCHASTIC WAVE-FUNCTION 
SIMULATIONS 



In this section wc perform Monte Carlo simulations of 
the generalized master equation for our model for both 
the weak coupling and strong coupling cases. The ter- 
minology, weak coupling and strong coupling are used in 
the same sense as described in Rcf. Details of the 
simulation algorithm can also be found in Ref. [l|. 



A. Weak Coupling 

It is clear to see that Eqs. ([15]) and pB]) are of the same 
form as Eq. ([1]) with the associations 



= = 0, R^^ = = 0, 



(24) 
(25) 



Here we have n = 2 and therefore consider two state 
vectors and \tp2)- The drift terms for the model 
from Eq. ^ are therefore given by 

Gi = -^(72'T+a_-72lk-|V'i)iri-7i|k+l^2)iri) 



^ [ 72 — 72C1 — 71 C2 

2 V -72C1 - 71C2 



with realizations 



\Mt)) 



-iGit 



^*IV'i 



-iGit 



l^l>ll 



(26) 



(27) 



and 



G2 = --(7i^-^+-72ik-iV'i)iri-7iik+iv^2)iri) 



^ / -72C1 - 71C2 

2 V 71 - 72C1 - 71 C2 



(28) 
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new TCL2 
new TCL2(t) 
Schrodinger 
MC 




with rate M'^ 



1000 



FIG. 2: (Color online) Comparison of the four methods with 
Ni ^ N2 = 200, Se = 0.31 and A = 0.001. 'new TCL2' 
and 'new TCL2(t)' correspond to Eq. (|36j and Eq. (|37)) 
respectively. The Monte Carlo simulation, 'MC, was done 
with time independent rates and the 'Schrodinger' gives the 
exact result. 



IV' 



72||a_|7^i)|p and 

<J+\iIJ2) 
lk+IV'2)||' 



(31) 



with rate AP = 7i||ct+|V'2)|P- 

The total waiting time distribution is 



F{t) = 1 - exp[- 



I^^a'IV',) 



exp[ 
exp[ 



El 

-7i||(T+|7/>2)|pr- 72||cr_|V'i)|pT] 
-7iC2r - 72Cit] (32) 



and depending on the current reaUzations, ci or C2 equals 
zero. It is easy to see that this process is rather simple, 
in that, beginning with the initial state |'0i(O)) = |e) and 
|V'2(0) = 0, the process simply jumps between = 
|e), |V2)=0and |^i) = 0, |^2> - Iff). 



B. Strong Coupling 



with realizations 



\Mt)) 



-tG2t 



1^2} 



-lG2t 



where ci = ||(T_|V'i)|P and C2 = ||cr+|V'2)|P- 
The two possible jumps arc 



IKI^OII 



(29) 



(30) 



In this case, Eqs. (fT5|) and are of the same form 
as Eq. Il]) with the associations 



(33) 
(34) 



^271 (T+, i?^-" V272Cr-. 

The drift terms and realizations are of the same form as 
for the weak coupling case, except here we need to take 
into consideration the time dependence in the waiting 
times. The total waiting time distribution is given by 




F{t) = 1 — exp[2 / c?ti/i(r — ti)(— 71C2T — 72C1T)] 
Jo 

"2(-l + cos(fcT) + Set Si(feT)) 



1 — exp 



Sett: 



(-71C2T - 72C1T) 



(35) 



r 



where Si(w) = ^^^dx. Once again, depending on the 
current realizations, ci or C2 equals zero. 



V. RESULTS 

In both cases we have considered the environment with 
Ni — N2 — 200 energy levels and the relaxation rates 7 ~ 
71 = 72- Se was chosen to be 0.31 so that for A — 0.001, 
0.013 and for A = 0.01, 4^ = 1.3. Note 



the ratio ^-^ 



71,2 

Se ~ 

that for the two cases considered, the relaxation rates 
differ by a factor 100. 

For the simulation of the new TCL2 with timc- 
indcpendcnt rates the waiting time distribution is 
F{ti^2) — cxp(— 72.1T1.2), which is just the exponential 



distribution, 
simulate 



For the initial condition pii{0) = 1, we 



Pii(i) 



-271 



(36) 



For the simulation of the new TCL2(t), the procedure 
is the same except that we need to include the time de- 
pendence in the waiting times. A Gaussian quadrature 
algorithm was used to evaluate the integral of /i(t — ti) 
and a polynomial interpolation algorithm was used to 
extract the waiting times, ti_2. With initial condition 
Pi 1(0) = 1, we simulate 



1 



(37) 
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FIG. 3: (Color online) Comparison of the four methods with 
Ni = N2 = 200, 5e = 0.31 and A = 0.01. 'new TCL2' 
and 'new TCL2(t)' correspond to Eq. (|36j and Eq. (|37)) 
respectively. The Monte Carlo simulation, 'MC, was done 
with time dependent rates and the 'Schrodinger' gives the 
exact result. 



where Ti^2{t) — 471,2 Jo dti J^^ dt2h{ti — 12). In both 
cases the Monte Carlo simulations were done with the 
initial state: |V'i(0)) = |e) and |V'2(0)) = 0. Also, in both 
cases, 5000 trajectories where used in the Monte Carlo 
simulations to recover the quantum master equation. 

We have also performed numerical solutions of the full 
Schrodinger equation corresponding to the Hamiltonian 
given in Eq. p^ . The initial state was taken to be 
|1) ® Ix), where the environmental state |x) was of the 
form 

(xl = (0r~^^~l0,di, ,d^J, (38) 

where di,...,dN-^ are Gaussian random variables with 
zero mean and variance equal to one. Ai?, the level dis- 
tance of the two-state system was taken to be unity. 



In Figs. [2] and [3] we compare the results of the four dif- 
ferent methods discussed in the paper, i.e., the new TCL2 
given by Eq. the new TCL(t) given by Eq. the 
numerical solution of the Schrodinger equation and the 
Monte Carlo simulations based on the unravelling of the 
master equation. For the weak coupling. Fig. [2] shows 
a good overlap of all four methods. For the strong cou- 
pling, as shown in Fig. O the Monte Carlo simulation 
results overlap almost completely with the new TCL2(t) 
method and also gives the correct stationary state and 
relaxation time when compared to the exact result ob- 
tained by solving the Schrodinger equation. 



VI. CONCLUSIONS 

In this paper, we have performed a stochastic unrav- 
elling of the generalized Lindblad master equation [2^ 
and applied it to a two-level system coupled to an envi- 
ronment consisting of two energy bands with 200 energy 
levels each. Our unravelling was applicable to both the 
weak coupling regime with time independent rates and 
the strong coupling regime with time dependent rates, 
for this model. Our Monte Carlo simulation results were 
found to be in good agreement with the second order 
time-convolutionless projection operator method results 
as obtained by the authors of Ref. [2^ . 
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